/* Angles in CIF are specified as an integer vector that has the proper
 * angle. These routines provide this funcionality.
 *
 * Copyright (c) 2006-2007 Xaveer Leijtens
 *
 * This file is part of libMask.
 *
 * libMask is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * libMask is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 * for more details.
 *
 * You should have received a copy of the GNU General Public License along
 * with libMask.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "mask.h"

extern double map_0_2pi(double);
extern void ToFract(double, int*, int*);

/* Return an integer vector (x,y) at specified angle a. Vector length is
 * not relevant.
 */
void angle_vector(double a, int *x, int *y)
{
  a = map_0_2pi(a);
  /* Special case: PI/2 and 3*PI/2 */
  if (fabs(a-M_PI/2) < 1e-10) {
    *x = 0; *y = 1; return;
  } else if (fabs(a-3*M_PI/2) < 1e-10) {
    *x = 0; *y = -1; return;
  }
  /* map to first quadrant: the ToFract function only works for positive
   * values.
   */
  if (a < M_PI/2) {
    ToFract(tan(a), x, y);
  } else if (a < M_PI) {
    ToFract(tan(M_PI-a), x, y);
    *x = -*x;
  } else if (a < 3*M_PI/2) {
    ToFract(tan(a-M_PI), x, y);
    *x = -*x;
    *y = -*y;
  } else {
    ToFract(tan(2*M_PI-a), x, y);
    *y = -*y;
  }
}

/* This function returns two integers. The division is an approximation of
 * the decimal number within a specified precision.
 *
 * It is needed for CIF files, where e.g. angles are specified as a
 * fraction of two integers.
 *
 * Modeled after http://mathforum.org/library/drmath/view/51886.html
 * (Doctor Peterson)
 */
void ToFract(double val, int *x, int *y)
{
  int lo[2], hi[2];
  /* 0=numerator, 1=denominator */
  double testlo;
  double testhi;
  double test;
  int intPart = (int)val;
  int i, count, num, denom;
  double Precision = 1e-12;

  val -= (double)intPart;
  lo[0] = 0; lo[1] = 1;
  hi[0] = 1; hi[1] = 1;
  for (i = 0; i < 100; ++i) {
    testlo = lo[1] * val - lo[0];
    testhi = hi[0] - hi[1] * val;
    if (testhi < Precision * hi[1]) break; /* high is answer */
    if (testlo < Precision * lo[1]) {	/* low is answer */
      hi[0] = lo[0]; hi[1] = lo[1];
      break;
    }
    if (i & 1) { /* odd step: add multiple of low to high */
      test = testhi / testlo;
      count = (int)test;
      num = (count + 1) * lo[0] + hi[0];
      denom = (count + 1) * lo[1] + hi[1];
/*      if ((num > 0x8000) || (denom > 0x10000)) break; */
      hi[0] = num - lo[0];
      hi[1] = denom - lo[1];
      lo[0] = num;
      lo[1] = denom;
    } else { /* even step: add multiple of high to low */
      test = testlo / testhi;
      count = (int)test;
      num = lo[0] + (count + 1) * hi[0];
      denom = lo[1] + (count + 1) * hi[1];
/*      if ((num > 0x10000) || (denom > 0x10000)) break; */
      lo[0] = num - hi[0];
      lo[1] = denom - hi[1];
      hi[0] = num;
      hi[1] = denom;
    }
  }
  *x = hi[1];
  *y = intPart * *x + hi[0];
}
